Specification Curve Analyses
Dani Cosme
2020-08-11
load packages
define palettes
pal_outcome = wesanderson::wes_palette("Zissou1", 3, "continuous")
pal_sca = c(pal_outcome[3], pal_outcome[3], "grey")
pal_condition = wesanderson::wes_palette("Zissou1", 16, "continuous")
pal_condition = c(pal_condition[1], pal_condition[3],
pal_condition[14], pal_condition[5],
pal_condition[15], pal_condition[16],
pal_condition[9], pal_condition[12])define SCA wrapper functions
plot functions
plot_sca = function(data, combined = TRUE, labels = c("A", "B"), title = FALSE, limits = NULL,
point_size = 1, point_alpha = 1, ci_alpha = .5, ci_size = .5,
text_size = 14, title_size = 6, n_breaks = 5,
choices = c("x", "y", "test", "stimulus", "contrast", "controls"),
remove_y = FALSE, remove_facet = FALSE) {
if (combined == TRUE) {
p1 = plot_curve(data, point_size = point_size, point_alpha = point_alpha,
ci_alpha = ci_alpha, ci_size = ci_size,
limits = limits) +
geom_hline(aes(yintercept = median(filter(data, y == "body fat percentage")$estimate), linetype = "body fat percentage"), color = pal_outcome[1], size = 1) +
geom_hline(aes(yintercept = median(filter(data, y == "BMI")$estimate), linetype = "BMI"), color = pal_outcome[2], size = 1) +
geom_hline(aes(yintercept = median(filter(data, y == "enactment")$estimate), linetype = "enactment"), color = pal_outcome[3], size = 1) +
scale_linetype_manual(name = "", values = c(2, 2, 2), guide = guide_legend(override.aes = list(color = c(pal_outcome[1], pal_outcome[2], pal_outcome[3])))) +
labs(x = "", y = "standarized\nregression coefficient\n") +
theme(legend.position = "top",
text = element_text(size = text_size))
if (title == TRUE) {
if (is.null(limits)) {
p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$x), fontface = 2, size = title_size,
x = 0.5*(1 + nrow(data)),
y = max(data$conf.high))
} else {
p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$x), fontface = 2, size = title_size,
x = 0.5*(1 + nrow(data)),
y = limits[2])
}
}
p2 = plot_choices(data, choices = choices,
ignore_vars = c("ROI", "neural signature"), rename_controls = "covariates") +
labs(x = "\nspecifications (ranked)") +
theme(strip.text.x = element_blank(),
text = element_text(size = text_size))
}
else {
p1 = plot_curve(data, point_size = point_size, point_alpha = point_alpha,
ci_alpha = ci_alpha, ci_size = ci_size) +
geom_hline(yintercept = 0, linetype = "solid", color = "black", size = .5) +
scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
labs(x = "", y = "standarized\nregression coefficient\n") +
theme(text = element_text(size = text_size))
if (title == TRUE) {
if (is.null(limits)) {
p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$y), fontface = 2, size = title_size,
x = 0.5*(1 + nrow(data)),
y = max(data$conf.high))
} else {
p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$y), fontface = 2, size = title_size,
x = 0.5*(1 + nrow(data)),
y = limits[2])
}
}
p2 = plot_choices(data, choices = choices,
ignore_vars = c("ROI", "neural signature", "not included"), rename_controls = "covariates") +
labs(x = "\nspecification number (ranked)") +
theme(strip.text.x = element_blank(),
text = element_text(size = text_size))
}
if (remove_y == TRUE) {
p1 = p1 + labs(y = "")
p2 = p2 + theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank()) +
labs(y = "")
}
if (remove_facet == TRUE) {
p2 = p2 + theme(strip.text.y = element_blank())
}
plot_specs(plot_a = p1,
plot_b = p2,
labels = labels,
rel_height = c(1, 2))
}
plot_sca_compare = function(data, pointrange = TRUE, labels = c("A", "B"),
rel_heights = c(.75, .25), rel_widths = c(.75, .25), b_limits = NULL,
title = FALSE, text_size = 14, title_size = 6, n_rows = 1, n_breaks = 3,
remove_x = FALSE, remove_y = FALSE, sig = NULL) {
# define median bootstrapping function
median_cl_boot = function(x, conf = 0.95, df = TRUE, ci = "low") {
lconf = (1 - conf)/2
uconf = 1 - lconf
require(boot)
bmedian = function(x, ind) median(x[ind])
bt = boot(x, bmedian, 1000)
bb = boot.ci(bt, type = "perc")
if (df == TRUE){
data.frame(y = median(x),
ymin = quantile(bt$t, lconf),
ymax = quantile(bt$t, uconf))
} else {
if (ci == "low") {
quantile(bt$t, lconf)
} else {
quantile(bt$t, uconf)
}
}
}
# merge and tidy for plotting
plot.data = data %>%
group_by(x) %>%
arrange(estimate) %>%
mutate(specification = row_number()) %>%
ungroup() %>%
unique()
# labels
labs = plot.data %>%
group_by(x) %>%
summarize(med = median(estimate),
low = median_cl_boot(estimate, df = FALSE, ci = "low"),
high = median_cl_boot(estimate, df = FALSE, ci = "high")) %>%
mutate(range = max(high) - min(low),
#estimate = ifelse(med > 0, high + (range / 10), low - (range / 10)),
estimate = ifelse(med > 0, high + .03, low - .03),
label = ifelse(x %in% sig, "*", ""))
# plot curves
if (pointrange == TRUE) {
a = plot.data %>%
ggplot(aes(specification, estimate, color = x)) +
geom_linerange(aes(ymin = conf.low, ymax = conf.high), size = .1) +
geom_point() +
geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 1) +
scale_color_manual(name = "", values = pal_condition) +
scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
labs(x = "\nspecification number (ranked)", y = "standarized\negression coefficient\n") +
theme_minimal() +
theme(strip.text = element_blank(),
axis.line = element_line("black", size = 0.5),
legend.position = c(.5, .1),
legend.direction = "horizontal",
panel.spacing = unit(0.75, "lines"),
axis.text = element_text(colour = "black"),
text = element_text(size = text_size))
if (title == TRUE) {
a = a + annotate("text", -Inf, Inf, label = unique(plot.data$y), fontface = 2, size = title_size,
x = 0.5*(min(plot.data$specification) + max(plot.data$specification)),
y = max(plot.data$conf.high))
}
} else {
a = plot.data %>%
ggplot(aes(specification, estimate, color = x)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 1) +
scale_color_manual(name = "", values = pal_condition) +
scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
labs(x = "\nspecification number (ranked)", y = "standarized\nregression coefficient\n") +
theme_minimal() +
theme(strip.text = element_blank(),
axis.line = element_line("black", size = 0.5),
legend.position = "none",
legend.direction = "horizontal",
panel.spacing = unit(0.75, "lines"),
axis.text = element_text(colour = "black"),
text = element_text(size = text_size))
if (title == TRUE) {
a = a + annotate("text", -Inf, Inf, label = unique(plot.data$y), fontface = 2, size = title_size,
x = 0.5*(min(plot.data$specification) + max(plot.data$specification)),
y = max(plot.data$estimate))
}
}
b = plot.data %>%
group_by(x) %>%
mutate(order = median(estimate)) %>%
ggplot(aes(reorder(x, order), estimate, fill = x)) +
stat_summary(fun = "median", geom = "bar") +
geom_errorbar(data = labs, aes(x = x, ymin = low, ymax = high), width = 0) +
geom_text(data = labs, aes(label = label, x = x, y = estimate), size = 6) +
scale_fill_manual(name = "", values = pal_condition) +
labs(x = "\n", y = "median effect\n") +
theme_minimal() +
theme(strip.text = element_blank(),
axis.line = element_line("black", size = 0.5),
legend.position = "none",
panel.spacing = unit(0.75, "lines"),
axis.text = element_text(colour = "black"),
text = element_text(size = text_size),
axis.text.x = element_text(angle = 45, hjust = 1))
if (!is.null(b_limits)) {
b = b +
scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
coord_cartesian(ylim = b_limits)
} else {
b = b +
scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks))
}
if (n_rows == 1) {
a = a + theme(legend.position = c(.5, .1))
b = b + coord_flip() +
labs(x = "\n", y = "\nmedian effect") +
theme(axis.text.x = element_text(angle = 0, hjust = 1),
axis.text.y = element_blank())
}
if (remove_x == TRUE) {
a = a + labs(x = "")
if (n_rows == 1) {
b = b + labs(y = "")
} else {
b = b + labs(x = "")
}
}
if (remove_y == TRUE) {
a = a + labs(y = "")
if (n_rows == 1) {
b = b + labs(x = "")
} else {
b = b + labs(y = "")
}
}
cowplot::plot_grid(a, b, labels = labels, rel_heights = rel_heights, rel_widths = rel_widths, nrow = n_rows)
}summary and bootstapping functions
Process overview
- use case resampling approach to generate a null dataset for each specification
- run all model specifications
- extract the dataset for each model specification
- force the null on each specification by subtracting the effect of the independent variable of interest (b estimate * x) from the dependent variable (y_value) for each observation in the dataset
- for each bootstrap iteration, sample with replacement from the null dataset for each specification and run the model specifications to generate a curve
- extract median estimate, % positive (& significant p < .05), % negative (& significant p < .05)
- repeat process 1000 times
run_boot = function(data, var, n_shuffles, rerun) {
if (rerun == FALSE & file.exists(sprintf("boot/boot_%s.RDS", var))) {
out = readRDS(sprintf("boot/boot_%s.RDS", var))
} else {
out = sca_boot(data = data, var = var, n_shuffles = n_shuffles)
saveRDS(out, sprintf("boot/boot_%s.RDS", var))
}
return(out)
}
sca_boot = function(data, var, n_shuffles, conf.level = 0.95) {
# set seed
set.seed(63)
# initialize data frame
df = data.frame()
# separate uniformity and association data
data1 = data %>%
filter(grepl("uniformity|ROI", test)) %>%
select(-test) %>%
spread(process, value_std)
data2 = data %>%
filter(grepl("association|ROI", test)) %>%
select(-test) %>%
spread(process, value_std)
# define control variables, subsets, and run scas
# define controls
if (length(var) > 1) {
controls = c(unique(data$process), "sex", "restraint")
} else {
controls = c(unique(data$process)[!unique(data$process) %in% var], "sex", "restraint")
}
# define subsets
subsets = list(stimulus = unique(data$stimulus),
contrast = unique(data$contrast))
# run enact_prop models separately
if (any(outcome == "enact_prop" & length(outcome) > 1)) {
# define controls
controls_enact = controls[controls != "sex"]
# run scas
results_data1_enact = run_specs(df = data1,
y = "enact_prop",
x = var,
controls = controls_enact,
model = c("lm"),
subsets = subsets,
keep.results = TRUE,
keep.formula = TRUE) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
results_data2_enact = run_specs(df = data2,
y = "enact_prop",
x = var,
controls = controls_enact,
model = c("lm"),
subsets = subsets,
keep.results = TRUE,
keep.formula = TRUE) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
}
if (any(outcome == "enact_prop") & length(outcome) == 1) {
# define controls
controls = controls[controls != "sex"]
# run scas
results_data1 = run_specs(df = data1,
y = outcome,
x = var,
controls = controls,
model = c("lm"),
subsets = subsets,
keep.results = TRUE,
keep.formula = TRUE) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
results_data2 = run_specs(df = data2,
y = outcome,
x = var,
controls = controls,
model = c("lm"),
subsets = subsets,
keep.results = TRUE,
keep.formula = TRUE) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
} else {
# run scas
results_data1 = run_specs(df = data1,
y = outcome[!grepl("enact_prop", outcome)],
x = var,
controls = controls,
model = c("lm"),
subsets = subsets,
keep.results = TRUE,
keep.formula = TRUE) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
results_data2 = run_specs(df = data2,
y = outcome[!grepl("enact_prop", outcome)],
x = var,
controls = controls,
model = c("lm"),
subsets = subsets,
keep.results = TRUE,
keep.formula = TRUE) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
}
# merge results dataframes
merged = bind_rows(results_data1, results_data2)
if (any(outcome == "enact_prop" & length(outcome) > 1)) {
merged = bind_rows(merged, results_data1_enact, results_data2_enact)
}
# create null dataset by subtracting the effect of x (estimate * x) from the dependent variable (y_value)
null_data = merged %>%
filter(!x == controls) %>%
mutate(stimulus = ifelse(grepl("snack", subsets), "snack",
ifelse(grepl("meal", subsets), "meal",
ifelse(grepl("dessert", subsets), "dessert",
ifelse(grepl("food", subsets), "food (average)", "all")))),
contrast = ifelse(grepl("nature", subsets), "nature",
ifelse(grepl("rest", subsets), "rest",
ifelse(grepl("average", subsets), "average", "all"))),
duplicated = ifelse(duplicated(estimate) & duplicated(std.error), 1, 0)) %>%
filter(!duplicated) %>%
unique() %>%
ungroup() %>%
filter(!stimulus == "all" & !contrast == "all") %>%
rownames_to_column() %>%
rename("model_number" = rowname) %>%
group_by(model_number) %>%
mutate(data = list(res[[1]][["model"]])) %>%
select(-res) %>%
unnest() %>%
group_by(model_number) %>%
mutate(obs_number = row_number()) %>%
ungroup() %>%
gather(outcome, y_value, bmi, fat, enact_prop) %>%
mutate(y_null = y_value - (estimate * !!as.name(var))) %>%
select(-y_value, -estimate) %>%
spread(outcome, y_null) %>%
select(-c(std.error, statistic, p.value, conf.low, conf.high, obs, obs_number, test))
for (i in 1:n_shuffles){
sampled = null_data %>%
group_by(model_number) %>%
mutate(bmi = ifelse(y == "bmi", sample(bmi, replace = TRUE), NA),
fat = ifelse(y == "fat", sample(fat, replace = TRUE), NA),
enact_prop = ifelse(y == "enact_prop", sample(enact_prop, replace = TRUE), NA)) %>%
group_by(model_number, x, y, model, formula, controls, random_effects, subsets, stimulus, contrast) %>%
nest()
results = sampled %>%
mutate(res = map2(.x = .data$data,
.y = .data$formula,
.f = ~ lm(.y, data = .x)),
coefs = map(.data$res,
broom::tidy,
conf.int = TRUE,
conf.level = conf.level),
obs = map(.data$res, nobs)) %>%
unnest(.data$coefs) %>%
unnest(.data$obs) %>%
filter(.data$term == .data$x) %>%
select(-.data$term, -.data$res, -.data$formula) %>%
mutate(pos = ifelse(estimate > 0, 1, 0),
neg = ifelse(estimate < 0, 1, 0),
sig = ifelse(p.value < .05, 1, 0),
pos_sig = ifelse(pos == 1 & sig == 1, 1, 0),
neg_sig = ifelse(neg == 1 & sig == 1, 1, 0)) %>%
group_by(x, y) %>%
summarize(median = median(estimate),
n = n(),
n_positive = sum(pos),
n_negative = sum(neg),
n_significant = sum(sig),
n_positive_sig = sum(pos_sig),
n_negative_sig = sum(neg_sig)) %>%
mutate(shuffle = i) %>%
select(shuffle, everything()) %>%
ungroup() %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "enact_prop", replacement = "enactment") %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "bmi", replacement = "BMI") %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "fat", replacement = "body fat percentage") %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "_", replacement = " ")
# join to data frame
df = rbind(df, as.data.frame(results))
rm(results)
}
return(df)
}
run_sca = function(data, var, outcome, boot = FALSE, n_shuffles = 1, rerun = FALSE) {
# set seed
set.seed(63)
# define median bootstrapping function
median_cl_boot = function(estimate, conf = 0.95) {
lconf = (1 - conf)/2
uconf = 1 - lconf
require(boot)
bmedian = function(estimate, ind) median(estimate[ind])
bt = boot(estimate, bmedian, 1000)
bb = boot.ci(bt, type = "perc")
data.frame(obs_median = median(estimate),
obs_conf_low = quantile(bt$t, lconf),
obs_conf_high = quantile(bt$t, uconf))
}
# separate uniformity and association data
data1 = data %>%
filter(grepl("uniformity|ROI", test)) %>%
select(-test) %>%
spread(process, value_std)
data2 = data %>%
filter(grepl("association|ROI", test)) %>%
select(-test) %>%
spread(process, value_std)
# define control variables, subsets, and run scas
# define controls
if (length(var) > 1) {
controls = c(unique(data$process), "sex", "restraint")
} else {
controls = c(unique(data$process)[!unique(data$process) %in% var], "sex", "restraint")
}
# define subsets
subsets = list(stimulus = unique(data$stimulus),
contrast = unique(data$contrast))
# run enact_prop models separately
if (any(outcome == "enact_prop") & length(outcome) > 1) {
# define controls
controls_enact = controls[controls != "sex"]
# run scas
results_data1_enact = run_specs(df = data1,
y = "enact_prop",
x = var,
controls = controls_enact,
model = c("lm"),
subsets = subsets) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
results_data2_enact = run_specs(df = data2,
y = "enact_prop",
x = var,
controls = controls_enact,
model = c("lm"),
subsets = subsets) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
}
if (any(outcome == "enact_prop") & length(outcome) == 1) {
# define controls
controls = controls[controls != "sex"]
# run scas
results_data1 = run_specs(df = data1,
y = outcome,
x = var,
controls = controls,
model = c("lm"),
subsets = subsets) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
results_data2 = run_specs(df = data2,
y = outcome,
x = var,
controls = controls,
model = c("lm"),
subsets = subsets) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
} else {
# run scas
results_data1 = run_specs(df = data1,
y = outcome[!grepl("enact_prop", outcome)],
x = var,
controls = controls,
model = c("lm"),
subsets = subsets) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
results_data2 = run_specs(df = data2,
y = outcome[!grepl("enact_prop", outcome)],
x = var,
controls = controls,
model = c("lm"),
subsets = subsets) %>%
mutate(test = ifelse(grepl("ROI", controls), "ROI",
ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI",
ifelse(controls == "craving_regulation_PEV", "neural signature",
ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
}
# merge results dataframes
merged = bind_rows(results_data1, results_data2)
if (any(outcome == "enact_prop" & length(outcome) > 1)) {
merged = bind_rows(merged, results_data1_enact, results_data2_enact)
}
# extract info & exclude duplicates
results = merged %>%
filter(!x == controls) %>%
mutate(stimulus = ifelse(grepl("snack", subsets), "snack",
ifelse(grepl("meal", subsets), "meal",
ifelse(grepl("dessert", subsets), "dessert",
ifelse(grepl("food", subsets), "food (average)", "all")))),
contrast = ifelse(grepl("nature", subsets), "nature",
ifelse(grepl("rest", subsets), "rest",
ifelse(grepl("average", subsets), "average", "all"))),
duplicated = ifelse(duplicated(estimate) & duplicated(std.error), 1, 0)) %>%
filter(!duplicated) %>%
unique() %>%
ungroup() %>%
filter(!stimulus == "all" & !contrast == "all") %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "enact_prop", replacement = "enactment") %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "bmi", replacement = "BMI") %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "fat", replacement = "body fat percentage") %>%
mutate_if(is.character, stringr::str_replace_all, pattern = "_", replacement = " ")
median_ci = results %>%
group_by(x, y) %>%
do({
median_cl_boot(.$estimate)
})
summary = results %>%
mutate(pos = ifelse(estimate > 0, 1, 0),
neg = ifelse(estimate < 0, 1, 0),
sig = ifelse(p.value < .05, 1, 0),
pos_sig = ifelse(pos == 1 & sig == 1, 1, 0),
neg_sig = ifelse(neg == 1 & sig == 1, 1, 0)) %>%
group_by(x, y) %>%
summarize(obs_n = n(),
obs_n_positive = sum(pos),
obs_n_negative = sum(neg),
obs_n_significant = sum(sig),
obs_n_positive_sig = sum(pos_sig),
obs_n_negative_sig = sum(neg_sig)) %>%
left_join(., median_ci, by = c("x", "y")) %>%
select(x, y, obs_median, obs_conf_low, obs_conf_high, everything())
if (boot == TRUE) {
boot_out = run_boot(data, var, n_shuffles, rerun)
return(list(results = results, summary = summary, boot = boot_out))
} else {
return(list(results = results, summary = summary))
}
}
boot_stats = function(data) {
data$boot %>%
mutate(y = ifelse(grepl("%", y), "body fat percentage", y)) %>%
left_join(., data$summary, c("x", "y")) %>%
mutate(extreme_median = ifelse(obs_median > 0 & median >= obs_median, 1,
ifelse(obs_median < 0 & median <= obs_median, 1, 0)),
extreme_positive = ifelse(n_positive >= obs_n_positive, 1, 0),
extreme_negative = ifelse(n_negative >= obs_n_negative, 1, 0),
extreme_positive_sig = ifelse(n_positive_sig >= obs_n_positive_sig, 1, 0),
extreme_negative_sig = ifelse(n_negative_sig >= obs_n_negative_sig, 1, 0)) %>%
group_by(x, y) %>%
summarize(n = n(),
extreme_median = sum(extreme_median),
extreme_positive = sum(extreme_positive),
extreme_negative = sum(extreme_negative),
extreme_positive_sig = sum(extreme_positive_sig),
extreme_negative_sig = sum(extreme_negative_sig)) %>%
gather(variable, n_extreme, contains("extreme")) %>%
mutate(p_value = n_extreme / n,
p_value = ifelse(p_value == 1.000, "1.000",
ifelse(p_value < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p_value))))) %>%
select(x, y, variable, p_value) %>%
spread(variable, p_value) %>%
left_join(., select(ungroup(data$summary), x, y, obs_median, obs_conf_low, obs_conf_high, obs_n_positive_sig, obs_n_negative_sig), by = c("x", "y")) %>%
mutate(`Mdn [95% CI]` = sprintf("%.2f [%.2f, %.2f]", obs_median, obs_conf_low, obs_conf_high)) %>%
select(y, `Mdn [95% CI]`, extreme_median, obs_n_positive_sig, extreme_positive_sig, obs_n_negative_sig, extreme_negative_sig) %>%
rename("Mdn p" = extreme_median,
"Positive N" = obs_n_positive_sig,
"Positive p" = extreme_positive_sig,
"Negative N" = obs_n_negative_sig,
"Negative p" = extreme_negative_sig) %>%
gather(var, val, -c(x, y)) %>%
unite(y, y, var, sep = " ") %>%
spread(y, val)
}load and tidy data
Read betas_dataset.RDS and dots_dataset.RDS saved from the xval_models.Rmd
- Ran uniformity and association test data separately because the models were breaking when trying to include test type as a covariate
- Code sex as -.5 (female), .5 (male)
source("load_data.R")
ind_diffs1 = ind_diffs %>%
mutate(bmi = as.numeric(scale(bmi)),
fat = as.numeric(scale(fat)),
restraint = as.numeric(scale(restraint)),
sex = ifelse(gender == 1, -.5,
ifelse(gender == 2, .5, gender))) %>%
select(-gender)
ema_enact1 = ema_enact %>%
mutate(enact_prop = as.numeric(scale(enact_prop)))
betas_std = readRDS("betas_dataset.RDS") %>%
filter(session == "all") %>%
group_by(subjectID, process, condition, control) %>%
mutate(value_std = mean(meanPE_std, na.rm = TRUE)) %>%
ungroup() %>%
select(subjectID, process, condition, control, value_std) %>%
unique() %>%
mutate(process = sprintf("%s_ROI", process)) %>%
filter(grepl("snack|meal|dessert|food", condition)) %>%
spread(control, value_std) %>%
mutate(average = (nature + rest) / 2) %>%
gather(control, value_std, nature, rest, average)
data_unif = readRDS("dots_dataset.RDS") %>%
filter(session == "all" & mask == "unmasked") %>%
filter((test == "uniformity" & !process == "craving_regulation") |
(test == "association" & process == "craving_regulation")) %>%
mutate(process = sprintf("%s_PEV", process)) %>%
rename("value_std" = dotProduct_std) %>%
select(subjectID, process, condition, control, value_std) %>%
unique() %>%
filter(grepl("snack|meal|dessert|food", condition)) %>%
spread(control, value_std) %>%
mutate(average = (nature + rest) / 2) %>%
gather(control, value_std, nature, rest, average) %>%
bind_rows(betas_std) %>%
left_join(., ind_diffs1, by = "subjectID") %>%
left_join(., ema_enact1, by = "subjectID") %>%
select(-c(sample, DBIC_ID, age)) %>% #gender
rename("contrast" = control,
"stimulus" = condition) %>%
mutate(stimulus = as.factor(stimulus),
contrast = as.factor(contrast)) %>%
spread(process, value_std) %>%
mutate(stimulus = gsub("food", "food (average)", stimulus))
data_assoc = readRDS("dots_dataset.RDS") %>%
filter(session == "all" & mask == "unmasked" & test == "association") %>%
mutate(process = sprintf("%s_PEV", process)) %>%
rename("value_std" = dotProduct_std) %>%
select(subjectID, process, condition, control, value_std) %>%
unique() %>%
filter(grepl("snack|meal|dessert|food", condition)) %>%
spread(control, value_std) %>%
mutate(average = (nature + rest) / 2) %>%
gather(control, value_std, nature, rest, average) %>%
bind_rows(betas_std) %>%
left_join(., ind_diffs1, by = "subjectID") %>%
left_join(., ema_enact1, by = "subjectID") %>%
select(-c(sample, DBIC_ID, age)) %>% #gender
rename("contrast" = control,
"stimulus" = condition) %>%
mutate(stimulus = as.factor(stimulus),
contrast = as.factor(contrast)) %>%
spread(process, value_std) %>%
mutate(stimulus = gsub("food", "food (average)", stimulus))
data_unif_gathered = data_unif %>%
gather(process, value_std, contains("ROI"), contains("PEV")) %>%
mutate(test = ifelse(grepl("PEV", process), "uniformity", "ROI"))
data_assoc_gathered = data_assoc %>%
gather(process, value_std, contains("ROI"), contains("PEV")) %>%
mutate(test = ifelse(grepl("PEV", process), "association", "ROI"))
data_combined = bind_rows(data_unif_gathered, data_assoc_gathered) %>%
unique()
head(data_combined)large coefficient models
specify modeling approach
- Also include the following subsets: stimulus type ( meal, snack, dessert, or food (average)) and contrast (nature, rest, or)
setup_specs(y = c("bmi", "fat", "enact_prop"),
x = c("reward_ROI", "reward_PEV",
"value_ROI", "value_PEV",
"cognitive_control_ROI", "cognitive_control_PEV",
"craving_PEV", "craving_regulation_PEV"),
controls = c("reward_ROI", "reward_PEV",
"value_ROI", "value_PEV",
"cognitive_control_ROI", "cognitive_control_PEV",
"craving_PEV", "craving_regulation_PEV"),
model = c("lm")) %>%
filter(!x == controls)define common variables
run SCAs
BMI
# define outcome
outcome = "bmi"
# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE)
bmi_output_plot = output$results %>%
mutate(`a priori` = ifelse(x == "value ROI" & controls == "reward ROI" &
stimulus == "food (average)" & contrast == "nature", "ROI model",
ifelse(x == "reward ROI" & controls == "value ROI" &
stimulus == "food (average)" & contrast == "nature", "ROI model",
ifelse(x == "craving regulation PEV" & controls == "no covariates" &
stimulus == "food (average)" & contrast == "nature", "PEV model", "not included"))))
# plot
plot_sca(data = bmi_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "controls", "a priori"))
% fat
# define outcome
outcome = "fat"
# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE, n_shuffles = n_shuffles, rerun = FALSE)
fat_output_plot = output$results %>%
mutate(`a priori` = ifelse(x == "reward ROI" & controls == "no covariates" &
stimulus == "food (average)" & contrast == "nature", "ROI model",
ifelse(x == "craving PEV" & controls == "no covariates" &
stimulus == "food (average)" & contrast == "nature" & test == "association", "PEV model", "not included")),
y = ifelse(grepl("%", y), "body fat percentage", y))
# plot
plot_sca(data = fat_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "controls", "a priori"))
enactment
# define outcome
outcome = "enact_prop"
# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE, n_shuffles = n_shuffles, rerun = FALSE)
enact_output_plot = output$results %>%
mutate(`a priori` = ifelse(x == "value ROI" & controls == "no covariates" &
stimulus == "food (average)" & contrast == "nature", "ROI model",
ifelse(x == "reward PEV" & controls == "no covariates" &
stimulus == "food (average)" & contrast == "nature" & test == "association", "PEV model", "not included")))
# plot
plot_sca(data = enact_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "controls", "a priori"))
combined plots
bmi_sca = plot_sca(data = bmi_output_plot, combined = FALSE, labels = c("A", "B"), title = TRUE,
point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
remove_facet = TRUE, n_breaks = 5,
choices = c("x", "test", "stimulus", "contrast", "a priori"))
fat_sca = plot_sca(data = fat_output_plot, combined = FALSE, labels = c("C", "D"), title = TRUE,
point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
remove_y = TRUE, remove_facet = TRUE, n_breaks = 5,
choices = c("x", "test", "stimulus", "contrast", "a priori"))
enact_sca = plot_sca(data = enact_output_plot, combined = FALSE, labels = c("E", "F"), title = TRUE,
point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
remove_y = TRUE, n_breaks = 5,
choices = c("x", "test", "stimulus", "contrast", "a priori"))
cowplot::plot_grid(bmi_sca, fat_sca, enact_sca, ncol = 3, rel_widths = c(.38, .31, .31))
bmi_compare = plot_sca_compare(data = bmi_output_plot, pointrange = FALSE, title = TRUE,
labels = c("A", "B"), rel_heights = c(.5, .5), text_size = 18,
n_rows = 2, n_breaks = 3,b_limits = c(-.25, .3),
sig = c("reward ROI",
"value ROI",
"cognitive control PEV",
"craving PEV", "craving regulation PEV"))
fat_compare = plot_sca_compare(data = fat_output_plot, pointrange = FALSE, title = TRUE,
labels = c("C", "D"), rel_heights = c(.5, .5), text_size = 18,
remove_y = TRUE, n_rows = 2, n_breaks = 3, b_limits = c(-.25, .3),
sig = c("value ROI",
"cognitive control ROI",
"craving PEV", "craving regulation PEV"))
enact_compare = plot_sca_compare(data = enact_output_plot, pointrange = FALSE, title = TRUE,
labels = c("E", "F"), rel_heights = c(.5, .5), text_size = 18,
remove_y = TRUE, n_rows = 2, n_breaks = 3, b_limits = c(-.25, .3),
sig = c("reward ROI", "reward PEV",
"value ROI", "value PEV",
"cognitive control ROI",
"craving PEV"))
cowplot::plot_grid(bmi_compare, fat_compare, enact_compare, ncol = 3)
neural covariate models
define common variables
run SCAs
reward ROI
# define variable
var = "reward_ROI"
# run SCA
reward_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
reward_ROI_output$summary 
reward PEV
# define variable
var = "reward_PEV"
# run SCA
reward_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
reward_PEV_output$summary 
value ROI
# define variable
var = "value_ROI"
# run SCA
value_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
value_ROI_output$summary 
value PEV
# define variable
var = "value_PEV"
# run SCA
value_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
value_PEV_output$summary 
cognitive control ROI
# define variable
var = "cognitive_control_ROI"
# run SCA
cognitive_control_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
cognitive_control_ROI_output$summary 
cognitive control PEV
# define variable
var = "cognitive_control_PEV"
# run SCA
cognitive_control_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
cognitive_control_PEV_output$summary 
craving PEV
# define variable
var = "craving_PEV"
# run SCA
craving_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
craving_PEV_output$summary 
plot combined SCAs
reward
limits = c(min(reward_PEV_output$results$conf.low), max(reward_PEV_output$results$conf.high))
a = plot_sca(data = mutate(reward_ROI_output$results,
controls = ifelse(controls == "reward PEV", "reward PEV / ROI", controls)),
combined = TRUE, title = TRUE, labels = c("A", "B"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_facet = TRUE, text_size = 18, limits = limits)
b = plot_sca(data = reward_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_y = TRUE, text_size = 18, limits = limits)
cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.54, .46))
value
limits = c(min(value_PEV_output$results$conf.low), max(value_PEV_output$results$conf.high))
a = plot_sca(data = mutate(value_ROI_output$results,
controls = ifelse(controls == "value PEV", "value PEV / ROI", controls)),
combined = TRUE, title = TRUE, labels = c("A", "B"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_facet = TRUE, text_size = 18, limits = limits)
b = plot_sca(data = value_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_y = TRUE, text_size = 18, limits = limits)
cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.53, .47))
cognitive control
limits = c(min(cognitive_control_ROI_output$results$conf.low), max(cognitive_control_PEV_output$results$conf.high))
a = plot_sca(data = mutate(cognitive_control_ROI_output$results,
controls = ifelse(controls == "cognitive control PEV", "cognitive control PEV / ROI", controls)),
combined = TRUE, title = TRUE, labels = c("A", "B"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_facet = TRUE, text_size = 18, limits = limits)
b = plot_sca(data = cognitive_control_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_y = TRUE, text_size = 18, limits = limits)
cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.57, .43))
craving & craving regulation
limits = c(min(craving_PEV_output$results$conf.low), max(craving_PEV_output$results$conf.high))
a = plot_sca(data = mutate(craving_PEV_output$results,
controls = ifelse(controls == "craving regulation PEV", "craving regulation / craving PEV", controls)),
combined = TRUE, title = TRUE, labels = c("A", "B"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_facet = TRUE, text_size = 18, limits = limits)
b = plot_sca(data = craving_regulation_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
choices = c("y", "test", "stimulus", "contrast", "controls"),
remove_y = TRUE, text_size = 18, limits = limits)
cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.57, .43))
combined table
bind_rows(reward_ROI, reward_PEV, craving_PEV, cognitive_control_ROI, cognitive_control_PEV, craving_regulation_PEV, value_ROI, value_PEV) %>%
kable(format = "pandoc", digits = 2)| x | BMI Mdn [95% CI] | BMI Mdn p | BMI Negative N | BMI Negative p | BMI Positive N | BMI Positive p | body fat percentage Mdn [95% CI] | body fat percentage Mdn p | body fat percentage Negative N | body fat percentage Negative p | body fat percentage Positive N | body fat percentage Positive p | enactment Mdn [95% CI] | enactment Mdn p | enactment Negative N | enactment Negative p | enactment Positive N | enactment Positive p |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| reward ROI | -0.13 [-0.14, -0.11] | < .001 | 22 | < .001 | 0 | 1.000 | -0.02 [-0.03, -0.01] | .012 | 0 | 1.000 | 0 | 1.000 | 0.19 [0.17, 0.21] | < .001 | 0 | 1.000 | 24 | < .001 |
| reward PEV | -0.00 [-0.01, 0.01] | .386 | 5 | .841 | 1 | 1.000 | -0.01 [-0.02, -0.00] | .045 | 4 | .930 | 0 | 1.000 | 0.21 [0.19, 0.22] | < .001 | 2 | .986 | 99 | < .001 |
| craving PEV | 0.04 [0.03, 0.05] | < .001 | 2 | .993 | 17 | < .001 | 0.04 [0.03, 0.05] | < .001 | 1 | .999 | 15 | .008 | 0.07 [0.05, 0.08] | < .001 | 7 | .486 | 23 | < .001 |
| cognitive control ROI | 0.02 [-0.00, 0.04] | .039 | 0 | 1.000 | 9 | .082 | 0.08 [0.07, 0.10] | < .001 | 0 | 1.000 | 12 | .007 | -0.15 [-0.19, -0.11] | < .001 | 34 | < .001 | 0 | 1.000 |
| cognitive control PEV | -0.07 [-0.07, -0.06] | < .001 | 22 | < .001 | 0 | 1.000 | -0.01 [-0.02, -0.00] | .002 | 2 | .991 | 0 | 1.000 | 0.00 [-0.01, 0.02] | .371 | 8 | .320 | 10 | .110 |
| craving regulation PEV | 0.09 [0.09, 0.11] | < .001 | 0 | 1.000 | 17 | < .001 | 0.06 [0.06, 0.08] | < .001 | 0 | 1.000 | 10 | .031 | -0.01 [-0.02, 0.00] | .246 | 1 | .993 | 0 | 1.000 |
| value ROI | 0.05 [0.03, 0.10] | < .001 | 0 | 1.000 | 56 | < .001 | 0.04 [0.03, 0.07] | < .001 | 0 | 1.000 | 35 | < .001 | 0.07 [0.05, 0.09] | < .001 | 0 | 1.000 | 19 | < .001 |
| value PEV | 0.02 [0.01, 0.04] | < .001 | 7 | .564 | 6 | .742 | 0.01 [0.00, 0.02] | .012 | 1 | 1.000 | 1 | 1.000 | 0.19 [0.17, 0.21] | < .001 | 0 | 1.000 | 81 | < .001 |